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A Particle-in-cell simulation model has been developed to study the physics of the Trav- 
eling Wave Direct Energy Converter (TWDEC) applied to the conversion of charged fusion 
products into electricity. In this model the availability of a beam of collimated fusion prod- 
ucts is assumed; the simulation is focused on the conversion of the beam kinetic energy 
into alternating current (AC) electric power. The model is electrostatic, as the electro- 
dynamics of the relatively slow ions can be treated in the quasistatic approximation. A 
two-dimensional, axisymmetric (radial-axial coordinates) geometry is considered. Ion beam 
particles are injected on one end and travel along the axis through ring-shaped electrodes 
with externally applied time-varying voltages, thus modulating the beam by forming a si- 
nusoidal pattern in the beam density. Further downstream, the modulated beam passes 
through another set of ring electrodes, now electrically floating. The modulated beam 
induces a time alternating potential difference between adjacent electrodes. Power can be 
drawn from the electrodes by connecting a resistive load. As energy is dissipated in the 
load, a corresponding drop in beam energy is measured. The simulation encapsulates the 
TWDEC process by reproducing the time-dependent transfer of energy and the particle 
deceleration due to the electric field phase time variations. 


Nomenclature 


a 

Acceleration, m/s 1 2 


Axial position, nr 

C 

Capacitance, F 

a 

Specific mass, kg/kW 

E 

Electric field vector, V/m 

Ah 

Spacing between grid nodes, nr 

TO 

Mass, kg 

At 

Time-step value, s 

P 

Power, W 

P 

Charge density, C/m 3 

q 

Ion bunch charge, C 

$ 

Electric potential, V 

Q 

Electrode charge, C 

Subscripts and superscripts 

r 

Radial position, nr 

i 

Radial grid node 

R 

Resistance, 

3 

Axial grid node 

V 

Velocity, m/s 

n 

Time-step 

V 

Velocity vector, m/s 

Physical Constants 

W 

Weight 

c 

Speed of light, 3 x 10 8 m/s 

X 

Position vector, m 

Co 

Permittivity of free space, 8.85 
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I. Introduction 


To change a crewed mission to Mars from a 3-year epic event to an annual (or sub-annual) expedition, 
new power and propulsion systems must be developed beyond the capabilities of both chemical rockets as well 
as nuclear thermal rocket designs. The key parameter to characterize the mission capability of a power and 
propulsion system is the specific mass, a, defined as the ratio of total mass to power output of the propulsion 
system (typically expressed in kg/kW). Under current Mars mission architectures (including nuclear thermal 
propulsion, solar electric propulsion, and solid core fission nuclear electric propulsion) a conjunction class 
mission to Mars has a manned mission duration of two to three years. 1 

Fusion concepts as well as advanced fission concepts have been proposed as the power system for electric 
propulsion systems. 2 Those concepts that release most of their energy in the form of neutrons require heavy 
shielding to protect the crew and systems from radiation, as well as a heat engine to convert the nuclear 
energy into electricity, necessitating large radiators to dispose of waste heat (figure 1 left). Aneutronic fusion 
and some types of fission, which release the majority of their energy in the form of charged particles, may 
utilize either heat energy conversion or the direct conversion of charged particles into electricity (such as 
in the TWDEC). Traditional heat energy conversion requires heavy radiators as well as additional power 
management and distribution hardware to condition the DC power output into AC suitable for the RF 
antenna in an ion engine. The TWDEC produces little waste heat, and can be physically tuned to match 
the frequency required for the ion engine RF antenna. By relaxing the mass requirements of radiation 
shielding, waste heat radiators, and DC to AC power management, a system using aneutronic fusion and a 
TWDEC may achieve a specific mass as low as 3 kg/kW (figure 1 center). 

The direct use of the momentum of fusion products for propulsion can further reduce specific mass. In 
this case the fusion products pass through the TWDEC en route to the direct plasma thruster. The TWDEC 
is used to decelerate the ions by only a small amount, to generate electricity for the spacecraft systems and 
maintain fusion core operation. The majority of the ion momentum is carried on to be used directly as 
thrust (figure 1 right). 
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Figure 1. A simplified comparison of the required component masses of three power and propulsion systems 
for a given jet power output. Box sizes roughly correspond to mass. 


II. The Physical TWDEC 

The principal mechanism behind the TWDEC is analagous to a linear particle accelerator working in 
reverse. Rather than putting energy into accelerating a beam, the TWDEC decelerates a beam and energy 
is extracted. The beam path may be considered in three sequential regions: the modulator, an intermediate 
region, and the decelerator, as displayed in figure 2 and outlined in the following sections. 

A. The modulator 

For the beam to become suitable for energy extraction it must be modulated from a homogeneous state into 
a modulated or “bunched” state. In the modulator section, a series of ring electrodes with time varying 
voltages create a wave of electric potential that travels with the beam, accelerating ions from the crests 
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Figure 2. Design for the NASA Johnson Space Center TWDEC test article. The ion beam travels left to 
right, through the modulator (A), intermediate region (B), and decelerator (C). 


towards the troughs, resulting in a roughly sinusoidal modulation of beam velocity while leaving unchanged 
the aggregate beam momentum. 

B. The intermediate region 

In the intermediate region, the perturbations in velocity space born in the modulator are given time to 
mature into density variations in the axial dimension, so that upon entering the decelerator section the ions 
are effectively “bunched” with a sinusoidal-like axial variation in density. 

C. The decelerator 

In the decelerator section the bunched ion beam passes through electrically floating electrode rings. The 
voltage of a floating electrode is determined by its surrounding environment; the proximity of a positively 
charged ion bunch will raise the floating electrode’s voltage. A traveling wave of sinusoidal beam density 
variation, then, will induce a sinusoidal voltage variation in a floating electrode. By adjusting the spacing 
of the electrodes, alternating voltage differences between electrodes arise. Power is produced by connecting 
these electrodes in a resistive circuit. 

A detailed explanation of bunch deceleration 

As the ion bunches induce time- varying voltages on electrodes, the ion bunch dynamics are in turn affected by 
the presence of these electrodes and any resistance in a circuit that might connect them. If passing through 
a single electrode, an ion bunch will experience an attraction to the electrode, due to the rearrangement of 
surface charge on the electrode that minimizes field energy. The ion bunch will accelerate axially towards the 
electrode as it approaches, and will decelerate accordingly after passing through, resulting in no net change 
in velocity and energy. * * 3. Similarly, if two electrodes are connected without resistance, no net deceleration 
occurs, as the charge transfer is immediate between electrodes, similar to an elastic collision. The two 
electrodes here act as one equipotential conductor (figure 3, middle). However, if a resistance is introduced 
between the electrodes, the transfer of charge is delayed, resulting in a net deceleration of the particle 
bunch. For the following steps, refer to the third plot in figure 3. When approaching the first electrode, the 
bunch is accelerated, but to a lesser degree than before, as the resistor delays the negative charge flow into 
that electrode (1). After passing through the first electrode, the bunch is decelerated (2). Upon reaching 
the midpoint of the electrode pair, electrode 1 remains positively charged while electrode 2 is negatively 
charged, due to the non-instantaneous nature of their discharging. At this point the deceleration continues, 
almost until the pair discharges (3), and the bunch undergoes a short acceleration in approaching the second 
electrode, reaching a local maximum in velocity (f). Upon passing through the second electrode the bunch 
undergoes deceleration due to an attractive force to the negatively charged electrode 2 (5). 

a This simple explanation neglects the expansion of the beam due to its space charge. In practice, radial beam expansion 

can be limited by an axial magnetic field. 
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Figure 3. Top: The axisymmetric 2-D domain with an ion bunch (here modeled as a semicircle which extrudes 
a sphere) passing through the centers of two ring electrodes (squares) Middle: A plot of bunch velocity vs. 
axial distance when the electrodes are connected by a wire with no resistance Bottom: The same plot with 
the electrodes connected by a 500 resistor, showing a net deceleration of the ion bunch. The relatively low 
magnitude of velocity loss is due to the pre-optimization stage of the model, as well as the lack of magnetic 
field to reduce bunch expansion. 


III. The Numerical Model 

The physics of the TWDEC is modeled with a particle-in-cell method that couples the Newtonian equa- 
tions of particle motion with the discrete Poisson equation to solve for the electric potential on cell nodes. 
The simulation exists in a 2-D axisymmetric radial-axial (r-z) domain, and it is assumed that motion and 
variance in the azimuthal dimension are negligible ( d/dd = 0). 

A. The particles 

The particles (ions) to be simulated are amalgamated into macroparticles to reduce computation times. The 
weight of each macroparticle, i.e. the number of particles each nracroparticle represents, is proportional to 
the inverse of the macroparticle’s initial distance from the axis upon generation. The nracroparticle weighting 
is only dependent on this initial radial distance; the weight of each macroparticle will remain constant over 
its lifetime in the simulation. The purpose of this non-uniform weighting is to have a uniform density of 
macroparticles in the 2-D plane. 
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B. Particle motion 


Particle motion is governed by the electrostatic equations of motion: 

dx civ q ^ 

- =v m = — e. 

at at m 

Discretizing the time variable into discrete time steps of value At, first order finite differencing is applied to 
the derivatives, with At kept small enough (in accordance with the CFL condition) to assume other variables 
stay constant over this time increment. 

-v-n+l ,,n+ 1/2 _ ,.n— 1/2 

X X _ v n+l/2 V V _ E n 

At At to 

Solving for the unknowns reveals the leapfrog method for particle transport. 

x” +1 = x" + v n+1/2 At v" +1/2 = v n -V 2 l E n At 

TO 

C. The matrix equation 

At each time step, given the charge distribution on the cell nodes, a system of equations is set up to 
simultaneously solve for the potential on the cell nodes, the surface charge distribution on the electrodes, 
and the current flow between electrodes. The borders of the floating electrodes exist on the cell nodes. 


1. Solving for the electric potential 

The potential is determined by the charge density in Poisson’s equation for electrostatics, 

V 2 <f> = -A 
£o 

In cylindrical coordinates with azimuthal symmetry ( d/d6 = 0), Poisson’s equation is 

< 9 2 $ 1 < 9 $ < 9 2 $ _ p 

dr 2 r dr dz 2 eo 

A first order finite difference approximation results in 

$;+ i,j ~ + $i-i ,j 1 ,j ®i,j + 1 - 2®i,j + $i,j - 1 = 

(A/i) 2 + iA/i 2A/i + (Ah) 2 ~ e 0 

simplifying to 

+ $j+iy(l + — ^) + + ^ij-i = -^-(A h)' 2 . 

For the special case when i = 0 (r = 0) it has been shown 3 that 

-6$ 0 y + 4$ij + $ 0j+ i + $ 0 ,i-i = - — (Ah) 2 . 

Co 

Thus for every discrete unknown $ in the domain there is a constraining equation. These equations can be 
put into matrix form and solved implicitly. 


2. Solving for the surface charge distribution on the decelerator electrodes 

The electrodes are approximated as perfect conductors, therefore the solver imposes the same potential on 
all nodes which reside on the electrode. Internal nodes are trivial and are neglected from the solver. The 
charge density on the surface node of the electrode is unknown and must be solved for, so these are added 
to the vector of unknowns. The conservation of charge on the surface provides an additional constraining 
equation, closing the system. 
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3. The circuit solution 

For electrodes connected in a circuit, two additional unknowns, the total charge on the each of the two 
electrodes, are introduced. Two additional constraining equations, then, are needed to close the system. As 
the potential difference builds between electrodes, current flows as in a discharging capacitor: 

dQ AV -t 

— = e«c . 

dt R 

Where C is the capacitance between the electrodes. This can be trivially approximated using finite differ- 
encing. 

, AV -At 

Q n = Q n 1 -| e RC At 

R 

The second constraining equation is the conservation of charge between the two electrodes: 
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Figure 4. The matrix system for the model, solving simultaneously for the electric potential at all nodes, the 
surface charge, and the total charge on each electrode 


The system can now be solved as an implicit matrix system. Solving is accomplished in MATLAB via the 
black box “backslash” operator. 

D. Particle-cell interpolation 

At each time step, values must be commuted from the particle locations to the cell nodes and vice versa. 
Quadratic interpolation 4 has been found to effectively avert the manifestation of non-physical furcation in 
beam expansion when lower order interpolation functions are used. The weighting scheme for quadratic 
interpolation is as follows. 

M<tt 

M 'M= kt-ik) 2 . M<" 

[O, otherwise 

The weighting schemes for each interpolation function are plotted in figure 5. A visual comparison of the 
effect of interpolation functions on an expanding beam can be seen in figure 6. 
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Figure 5. A graphical comparison of the different interpolation weighting methods. 



NEAREST NEIGHBOR LINEAR QUADRATIC 


Figure 6. A qualitative comparison of the effects of different interpolation methods on an expanding beam. 
A 100 X 50 grid is used, with approximately 33,000 particles. The chosen interpolation method in each case 
is used for both particle-to-grid as well as grid-to-particle interpolation. 


IV. The Quasistatic Approximation 


The justification for a quasistatic treatment of ion dynamics is dependent upon proof that the electro- 
magnetic energy radiation from a decelerating ion or ion bunch is negligible in comparison to the magnitude 
of direct energy conversion. For a test condition, an ion bunch losing 10% of its velocity through a pair of 
decelerator rings is considered. 15 Given a p-B 11 a-particle fusion product with a velocity of 10' m/s and a 
spacing between electrodes of 5 cm, an ion will experience a deceleration on the order of 10 11 m/s 2 . 


1 (An) 2 


1 (10 6 m/s) 2 
2 ~ 


= 10 u m/s 2 


2 Ax 2 .05 nr 

The power of this wave (assuming an ion bunch charge of 10 -9 C) can then be calculated using the Larmor 
formula 

1 q 2 a 2 1 (10 -9 C) 2 (10 11 m/s 2 ) 2 


P = 


10“ n W 


67T£o C 3 67T£o c 3 

The deceleration will take place over a time of approximately 10 -8 s, resulting in a radiated energy on the 
order of 10“ 19 J. Simulation results have shown that a similarly charged bunch will lose on the order of 
1CU 8 J when undergoing deceleration between two electrodes. The radiated energy of a decelerating bunch 
is many magnitudes lower than the electrostatic energy transfer, therefore the quasistatic approximation is 
valid for the model. 


V. Computational Acceleration with GPU Parallelization 

Computation times are reduced through partial execution on a high performance graphics processing 
unit (GPU). Parallelization is achieved through Jacket in MATLAB. Jacket provides black box operation 
for running iteration-independent loops and vectorized operations in parallel on the GPU, resulting in sig- 
nificant speed increases over high-end CPU processing. While some operations in the PIC code benefit 
prodigiously from GPU parallelization, other operations are unable to be executed on the GPU under the 
current capability of the software. In particular, the large matrix inversion operation involved in solving the 
Poisson equation, as well as a portion of the charge distribution interpolation, are performed on the CPU. 
A breakdown of these computation times is provided in figure 7. 

b It is not expected for a single pair of decelerator rings to ever achieve this magnitude of deceleration 
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Figure 7. A comparison of the amount of time required to execute a single time step in a simulation on a 
1000 x 375 grid with 1,000,000 particles. The CPU is a 12 core (2x6) 2.8 GHz Xeon X5660 and the GPU is 
an NVIDIA Tesla C2070. 


VI. Model Capabilities 

The PIC model conceived of and developed by the authors is currently capable of modeling the modulation 
and deceleration of an ion beam. The model has been partially validated by noting the conservation of 
energy between the beam kinetic energy, the electric field energy, and the resistive load energy. Further 
validation was provided by comparison of various model sub-physics with replications in COMSOL. The 
model is suitable for investigating the relation between energy conversion efficiency and electrode shape and 
positioning. Although past computational 5 and experimental 6 studies on the TWDEC have investigated the 
achievable beam deceleration efficiency, no studies have investigated computationally or experimentally the 
energy conversion mechanism at play in these devices. 
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Figure 8. Numerical simulation of the TWDEC with heuristic parameters. Dots denote macroparticle loca- 
tions, boxes denote electrodes, and lines of constant potential are plotted in color. The beam travels from left 
to right. Four modulator rings apply a traveling wave in potential with a wave velocity equal to the beam 
velocity. The downstream decelerator ring pairs experience alternating potentials. 


VII. Future Work 

Future work will include improvements to the model design, as well as optimization routines to direct 
the construction of an experimental TWDEC test article at NASA’s Johnson Space Center. 

Model improvement 

To improve the physics of the model, simulation of the collision of ions with background gas will be introduced, 
which will allow to model remain accurate for a laboratory vacuum. Because an integral part of TWDEC 
operation is the magnetic field to mitigate beam expansion, the motion of the gyrocenters of macroparticles 
in a magnetic field will be modeled. A proposed method of neutralization of the ion beam center with an 
electron beam will also be considered. 
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TWDEC optimization 

The primary motivation for this work is the optimization of the TWDEC specific mass. Of particular interest 
is the way in which electrode design and spacing affects energy transfer efficiency, as well as what efficiencies 
are achievable for a given beam energy and density. The model will also be used for future vehicle studies 
for which the TWDEC may be considered. 
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